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Abstract 



In this paper we propose a decentralized sensor network scheme capable to reach a globally optimum maximum 
likelihood (ML) estimate through self-synchronization of nonlinearly coupled dynamical systems. Each node of 
' the network is composed of a sensor and a first-order dynamical system initialized with the local measurements. 
, Nearby nodes interact with each other exchanging their state value and the final estimate is associated to the state 



derivative of each dynamical system. We derive the conditions on the coupling mechanism guaranteeing that, if the 
network observes one common phenomenon, each node converges to the globally optimal ML estimate. We prove 



CZ2 

J-J ' that the synchronized state is globally asymptotically stable if the coupling strength exceeds a given threshold. 



Acting on a single parameter, the coupling strength, we show how, in the case of nonlinear coupling, the network 



5h ' behavior can switch from a global consensus system to a spatial clustering system. Finally, we show the effect of 
the network topology on the scalability properties of the network and we validate our theoretical findings with 
simulation results. 



1 Introduction 

Sensor networks are receiving a significant attention because of their many potential civilian and military 

applications (see, e.g., [U O |3]). The single major challenge is perhaps how to conjugate the relative 

unreliability of the single node, due to its limited complexity and energy availability, with the high 

reliability required to the whole network. Most research works aim then at making the best use of the 

*This work has been partially funded by ARL, Contract N62558-05-P-0458, and by the WINSOC project, a Specific 
Targeted Research Project (Contract Number 0033914) co-funded by the INFSO DC of the European Commission within 
the RTD activities of the Thematic Priority Information Society Technologies. 

1 



available resources. Many works concentrate on how to adapt the protocol stacks derived in decades of 
research in communication networks to the sensor network scenario. An alternative approach consists, 
instead, in recognizing that a sensor network is intrinsically different from a communication network, 
thus implying that the design of a sensor network should reflect its specificities. Among the features 
distinguishing a sensor network from a communication network, we may mention its data-centric and 
event-driven nature. In general, a sensor network can be seen as a sort of distributed computer that, 
on the basis of the measurements, let us say xi,X2, ■ ■ ■ ,xn, gathered by sensors, it has to take a 
decision about the observed phenomenon by computing a function f{xi,X2, rcAr) of the measurements. 
Typically, this function has properties, depending on the application, that, if properly exploited, can 
suggest efficient ways to design a sensor network. For example, Giridhar and Kumar recently proved that, 
if /(xi,X2, . . . ,X]sr) is invariant to any permutation of the observed variables (like in the computation 
of the average, for example, or the maximum, etc.), it is possible to improve the scalability properties 
of the network, using some kind of in-network processing [4J. Interestingly, this symmetry property is 
not at all artificial, as it reflects the data-centric nature of the network and it holds true in a variety 
of applications. Some works exploit the data-centric property to devise innovative schemes, like, for 
example, the type-based multiple access (TBMA) system [6j. TBMA is perfectly scalable, but it requires 
a high coherence of the channels from the sensors to the sink node. 

Quite recently, several authors have proposed an alternative approach that allows each node to per- 
form in-network processing, so as to reduce the burden of the fusion center pj. Other works go even 
further by proposing strategies where the global decision, or estimation, is obtained using a totally dis- 
tributed approach, with no need for a fusion center, at least in the case where the whole network observes 
a common event (Tj (H [U [TOl [TTl [T2l [T9] . A strategy that has received significant attention in the last 
few years is the so called average consensus protocol. The basic idea is that, if the network is connected, 
i.e., there is a path, possibly composed of multiple hops, between any pair of nodes, local exchange of 
information among nearby sensors is sufficient to reach a global consensus on the average of the observed 
values, without requiring any control node. Global consensus can be reached through linear coupling, as in 
[19], [9], or through nonlinear coupling, as in [13], [H], [3l]. Global consensus can also be used to track a 
common time- varying phenomenon, as in [17], [16]. An important synergism to this approach comes also 
from the algorithms developed for the coordination of groups of mobile autonomous agents through local 
transmissions [2lj. An alternative approach to achieve a consensus was proposed in [22 | 1231 [M] . where 
consensus was seen as a result of self-synchronization of a population of pulse-coupled oscillators, each one 
initialized with the sensor local estimates or decisions. The principle ensuring the self-synchronization 
capability of the system proposed in [22], [23] relied on a theorem, proved by Mirollo and Strogatz in 
|27j . that required the full network connectivity, i.e., the property that each node has a direct link to 
each other node. This assumption was later removed by Lucarelli and Wang in [24J, who proved that 
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local coupling among the nodes is sufficient, provided that the whole network is connected. The pulse 
coupling mechanism is indeed appealing from the implementation point of view, but, especially for large 
scale networks, it may suffer from ambiguity problems, as the information bearing time shift may become 
indistinguishable from the propagation delay. The idea of achieving global estimates or decisions exploit- 
ing local coupling among dynamical systems, initialized with local measurements, was then proposed in 
|12j , [25] , [26] . Average consensus through mutual coupling of first order dynamical systems was used in 
|12j . |26| to derive the globally optimal maximum likelihood estimation: In [12] , each system is initialized 
with the local measurement, the coupling is linear and the consensus amounts to requiring all dynamical 
systems to reach the same value of the state; in [26], each system is randomly initialized, the coupling 
is nonlinear (it subsumes linear coupling as a particular case) and the consensus refers to the situation 
where the derivatives of the states (rather than the states) converge, asymptotically, to a common value. 

This paper builds on the initial idea of [26] and its main contributions are the following: i) we derive 
the conditions guaranteeing that a globally optimum maximum likelihood estimator can be obtained 
through local nonlinear coupling of first order dynamical systems; ii) we show that nonlinear coupling 
offers a variety of behaviors, to be used to find out the best implementation of the radio transceivers or 
to allow the network to work as a global estimator or as a spatial clustering mechanism; iii) we show that 
convergence on the state derivative (rater than the state) improves the resilience against additive noise 
with respect to common average consensus techniques. 

The paper is organized as follows. In Section [2] we describe the coupling mechanism. In Section [3] we 
show how to design the coupling mechanism and the local processing in order to make the equilibrium 
achievable by each node to coincide with the globally optimum maximum likelihood estimate. In SectionH] 
we derive the conditions guaranteeing that the equilibrium is unique and asymptotically stable. Finally, 
in Section [5] we report numerical results validating our theoretical findings and showing the network 
behavior both as a global estimator or as a spatial clustering system. 

2 Coupling mechanism 

The proposed sensor network is composed of N nodes, each composed of four basic components: i) a 
transducer that senses the physical parameter of interest (e.g., temperature, concentration of contami- 
nants, radiation, etc.); ii) a local detector or estimator that processes the measurements taken by the 
node; iii) a dynamical system whose state evolves according to a first-order differential equation, whose 
parameters depend on the local estimate and on the states of nearby nodes; iv) a radio interface that 
transmits the state of the dynamical system and receives the state transmitted by nearby nodes, thus 
ensuring the interaction among nearby nodes. 



2.1 Scalar observations 

When each sensor measures a single physical parameter, the dynamical system present in node i evolves 
according to the following equation 

K ^ 

^S=i i = l,...,iV, (1) 

6'i(0) = ^io, 

where 

1. 9i{t) is the state function of the z-th sensor, initialized, at f = 0, as any random number ^i(O) = 6iQ] 

2. cjj is a function g{xi) of the observation Xi taken from node i; 

3. /(•) is a nonlinear, odd function that takes into account the coupling among the sensorj^; 

4. K \s a. positive control loop gain measuring the coupling strength; 

5. Cj is a positive coefficient that quantifies the attitude of the i-th sensor to adapt its state as a 
function of the signals received from the other nodes; 

6. the coefficients Ojj take into account the local coupling among the systems: if nodes i and j are 
coupled to each other, aij ^ 0, otherwise Ojj = 0; we assume that the nonzero coefficients aij are 
positive and respect the symmetry condition a^j = aji; 

7. Vi{t) is additive noise. 



The mode 



([T]) coincides with the so called Kuramoto model [28], when /(x) = sin(x) and a^j = 
Given the model ([1]), the running decision, or estimate, of each sensor is associated 
to the derivative of the state function Oi{t). Global consensus, in this paper, means that all nodes 
end up evolving with the same state derivative. This choice is different from common average consensus 
techniques, where the consensus refers to the state value. We will show, in Section 5, that this apparently 
slight difference brings important consequences in the presence of additive coupling noise. Furthermore, 
in our case, the states of different nodes are let free to converge to functions differing by a constant term. 
This extra degree of freedom, with respect to the techniques converging on the state, might be exploited 
for different scopes than synchronization, like, for example, spatial pattern recognition, as in |30j . 
A possible schematic implementation of our protocol is reported in Fig. [TJ On the right side, there is 
a transducer that measures a physical quantity Xj and produces a parameter = g{xi). On the left 



^We assume, w.l.o.g., that f{x) is normalized so that /'(O) = 1, where /'(x) := df{x)/dx, as different values of /'(O) can 
always be included in K. 

^As will be clarified in Section 4, our function f{x) has to be a monotonically increasing function, to guarantee the 
achievement of the global optimal ML estimate. Hence, Kuramoto model is mentioned here only for similarity reasons, but 
our main findings do not include Kuramoto model as a particular case. 
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Figure 1: General coupling mechanism. 



side, there is a single transmit/receive antenna and a circulator used to switch between transmission and 
reception. The transmitted signal is a waveform p{t;6i{t)) that depends on the local state 6i{t). The 
received signal is a linear combination of the signals transmitted by nearby nodes, i.e. 

N 



where the coefficients depends on the propagation and radio interface. The received signal is then 
mixed with a local waveform q{t; 6i{t)), generated as a function of 9i{t), in order to produce the signal 

N 

y{t) = ^a,,f[9,{t)-9,{t)]. 
i=i 

It is easy to check that the input of the waveform generator coincides with ([T]) . 

The scheme of Fig. [T]is rather general and it incorporates alternative implementations, like pulse coupled 
systems, for example, or phase- locked loops, depending on the choice of the waveforms p{t) and q{t). The 
half-duplex feature of the scheme depicted in Fig. [T] implies that, if two nodes i and j transmit the same 
waveform, i.e. 9i{t) = 9j{t), they do not listen to each other. This feature is consistent, in mathematical 
terms, with the odd property of the function f{x), as will be described in Section 4. 

We assume, initially, that the additive noise is negligible, i.e. Vi{t) = 0. In Section El we will show 
the effect of noise on the system performance. 

To make explicit the network connectivity properties, it is useful to rewrite ([1]) introducing the graph 
incidence matrix B, defined as follows. Given an oriented graph composed by N vertices and E 



edges, B is the N x E matrix such that [B]ij = 1 if the edge j is incoming to vertex i, [B] 



-1 if the 



edge j is outcoming from vertex i, and otherwise. Given the x 1 vector Iat, composed of all ones, it 
is easy to check that the incidence matrix satisfies the following property: 



llB = 0%. 



(2) 



Given B, the symmetric N x N matrix L defined as -L = BB^ , is called the Laplacian of and it is 
independent of the graph orientation. If we associate a positive number Wi to each edge and we build 



The orientation of a grapli consists in the assignment of a direction to each edge. 



the diagonal matrix = diag(w), with w = [wi, ■ • • ,we]'^, we may introduce the so called weighted 
Laplacian, which is written as — BD^B"^. The Laplacian (as well as the weighted Laplacian) 
has several important properties, among which [3T]: 1) L (or L^) is always positive semi-definite, i.e. 
with the smallest eigenvalue always equal to zero; 2) the algebraic multiplicity of the null eigenvalue is 
equal to the number ric of connected components of the graph; if the graph is connected, Uc = 1 and 
rank(L) = rank(Lw) = A^ — 1, i.e. L (or Lw) has a unique zero eigenvalue and the eigenvector associated 
to the null eigenvalue is the vector Iat. The second smallest eigenvalue A2(L) (or A2(Lw)) is known as 
the graph algebraic connectivity, and it provides a measure of connectivity [32]. 

Using the above notation, La — BDaB"^ will denote the weighted Laplacian associated to the graph 
describing our network ([1]), including the positive coefficients {aij}. Furthermore, dmax — niaxj Ylf=i ^ij 
will denote the maximum degree of the (weighted) graph. 
Using the incidence matrix we can rewrite ([I]) in compact form as 

9{t) =u-KD-^BDAf [B^eit)] , 

("Jj 

e{o) = Oo, 

where 6{t) = [9i{t), ■ ■ ■ ,9N{t)]'^; 9q = [9iq, ■ ■ ■ ,9]\fo]'^; = diag {ci, . . . , cat}; Da is an ii^ x E' diagonal 
matrix, whose diagonal entries are all the weights Ojj, indexed from 1 to E] the symbol f{x) has to be 
intended as the vector whose k-ih. component is /(x^). 

2.2 Vector Observation 

When each sensor measures more, let us say L, physical parameters like, e.g., temperature, pressure, etc., 

J3 



the coupling mechanism ([T|) generalizes according to the following expressio: 

N 



e,{t) =u, + KQT^ V aij f [Ojit) - Oiit)] 



j=i i = l,...,N, (4) 

^j(O) = 9io, 

where Oi{t) is the L-size vector state of the i-th node, that is initialized as a random vector OiQ] Ui is the 
L-size vector, function of the L measurements taken by node i; Q, is an L x L non-singular matrix that 
depends on the observation model. In Section 3, we show how to choose the vectors cjj and the matrices 
Qj to guarantee the convergence of ([4]) to the global optimal maximum likelihood (ML) estimate. Also in 
this case, we can rewrite dH) in compact form using the graph incidence matrix. Introducing the vectors 
0{t) = [0f (t), . . . 6'^{t)]^ and u = [o;^, . . . , and the matrix Dq = diag(Q;^, . . . , Q^), the system 

(jl]) becomes 

e{t) =u-KT>^ P^(/i ® ED a) f [{II S^) P 0{t)] , 

(5) 

0(0) = Oo, 



''We assume that the coupling coefficients atj are the same for all estimated parameters. This assumption is justified by 
the fact that aij depends on the coverage radius of each transmitter and not on the measurements. 



where (8" denotes the Kronecker product, 0(0) = [6^(0), . . . Ojf{0)]'^ and P is an LN x LN permutation 
matrix defined as 

1, if j = (((imodiV) - 1)L + 
0; otherwise, 



I 

N 



mod( NL), 

(6) 



so that, in ([5]), each vector x = [x^, . . . ,x^]^, with Xj = [x\^\ . . . ,x\^^]'^, is mapped into a new vector 



X = Px, partitioned as x = [x-,^ , . . . , xj^] , with Xj = [x\ 



2.3 Network Self-Synchronization 



Differently from [12], where the global consensus was intended to be the situation where all dynami- 
cal systems reach the same state, we adopt an alternative consensus strategy. We define the network 
synchronization with respect to the state derivative, as follows: 

Definition 1 The overall population of dynamical systems ( or ^ ) is said to synchronize if there 
exists a solution 0*{t) of ([ip (or ^) such that all the state derivatives converge, asymptotically, to 0* (t), 
i.e. 

Ihn WOiit) - e\t)\\ = 0, i = l,2,...,iV, (7) 

ti— >oo 

where \\ ■ \\ denotes some vector norm. The system is globally asymptotically stable if the system syn- 
chronizes, in the sense specified before, for any set of initial conditions Oi{0). 

According to Definition [H if there exists a synchronized state that is globally asymptotically stable, then 
it must necessarily be unique. Interestingly, if the synchronized state exists, it can be computed in closed 
form, without explicitly solving the system of differential equations ([1]) and In fact, exploiting the 
oddness property of f{x), left-multiplying 1^ by the row vector = l^Dc, we obtain 

c^e{t) = - K I^BDa/ (B^0) = Ju, (8) 

where in the second equality of ([8]), we have used Hence, if system ([3]) synchronizes (according to 
Definition [1]) , the common value of {t) must be constant and equal to 



W=-^ = ^ = ^^^^- (9) 



Similarly, in the vector case, left-multiplying ([5]) by the matrix (l^ (gi 1^) Dq, we obtain 

N N 

^Qieiit) = ^QiUi, (10) 

i=l i=l 

where we used the following chain of equalities (1^ 1l)'P'^ (1^ ® BDa) = (II ® 1^)(Il ® BDa) = 
Il^I^BDa = 0l xLEi and the property ([2]). Hence, if the system synchronizes, in the sense of Definition 



1, the synchronized state must necessarily be 

3 Reaching global ML estimate through self-synchronization 

The basic idea of this paper is that, when the whole network observes a common phenomenon, the 
self-synchronization process forms the basic mechanism for reaching the globally optimal ML estimate 
through local exchange of the state functions, without sending the observations to any fusion center. In 
particular, let us consider the scalar observation 

Xi = biC + Wi, (12) 

where ^ is the common unknown parameter to be estimated and Wi, i = 1,...,N are a set of i.i.d. 
Gaussian random variables with zero mean and variances {af,...,a'j^). Initializing each node with 
uJi = Xi/bi and setting, in ([3]), Cj = 6? /cr?, the network synchronized state Q becomes: 

= ^M. = (13) 

This value coincides with the globally optimal ML estimate [33]. The equilibrium ()13p shows that the 
most reliable nodes (i.e., the ones with the smallest Oi) are the most influent nodes in driving the whole 
system towards the common decision. What is important to stress is that this happens without any node 
knowing which are the nodes with the best SNR. 

In the linear vector case, each node observes the vector 

■x.i = Aii^Wi, (14) 

where Xi is the M x 1 observation vector, ^ is the L x 1 unknown common parameter vector, Aj is 
the M X L mixing matrix, and Wi is the observation noise vector, modeled as a Gaussian vector with 
zero mean and covariance matrix Cj. We assume that the noise vectors affecting different sensors are 
statistically independent of each other (however, the noise vector present in each sensor may be colored). 
We consider the case where the single sensor must be able, in principle, to recover the parameter vector ^ 
from its own observation. This requires that M > L and that Aj is full column rank. In this case, setting, 
in ([5]), Qi = AfC^^Ai and Ui = [A^ C^"^ Aj)"^ AfC^^xi, the synchronized state ([U]) becomes 

^1 = 1ml = ^^Cr^^}j (j2 AfCr^x}j , (15) 
which coincides with the globally optimal ML estimate [33j. 



It is important to emphasize here that a virtually optimal fusion center in this case would need to 
know not only the observation vectors Xi, but also all mixing matrices Aj and the noise covariance 
matrices Cj. Conversely, following the proposed approach, if the network converges, each node tends 
to the optimal ML estimate without sending all these data to any sink node, but simply exchanging 
the state vectors 6i{t) with nearby nodes. The penalty for having this advantage is that the solution is 
reached through an iterative procedure that consumes time and energy. However, this energy is spent 
only for local transmissions. The crucial point, as far as the global energy consumption is concerned, is 
then the convergence time. In the next section, we will give an approximate formula for such a value. 

The proposed system has indeed a broader applicability than just global ML estimation for a linear 
observation model. From ([9|), since uJi is a function gi{xi) of the local observation, our approach allows 
us to compute any function of the collected data expressible in the form 



with positive coefficients Cj, with a totally distributed mechanism. Of course, this class of functions is not 
the most general one. Nevertheless, the class of functions in ()16p contains not only the linear ML case, 
but it comprises many cases of practical interest (like, for example, computation of the sufficient statistics 
in detection of Gaussian processes in Gaussian noise, computation of maximum, minimum, histograms, 
geometric mean of the observed values, etc.), as it can be checked by choosing the functions g{x),h{x) 
and the coefficients Cj appropriately. 

4 Global asymptotic stability of the synchronized state 

Given the dynamical system ([I]) or (jH), the natural questions to ask are: i) Does the synchronized state 
exist? ii) If it exists, does the system synchronize for any set of initial conditions? In this section, we 
provide an answer to these questions. We focus, initially, on the scalar system ([T]) and provide necessary 
and sufficient conditions for the existence of a globally asymptotically synchronized state, according to 
Definition 1. Then, we generalize the result to the vector system (jl]). 

Theorem 1 Given the system ([7]j, assume that the following conditions are satisfied: 
al The graph associated to the network is connected; 

a2 The nonlinear function /(•) : M i-^ M zs a continuously difjerentiable, odd, increasing function; 
a3 The nonzero coefficients Uij and the coefficients Ci are positive. 

Then, there exist two unique critical values of K, denoted by Kl and Kjj, with < < Kjj, such 
that the synchronized state exists for all K > Kjj, and it does not for all K < K^. Furthermore, if it 




(16) 



n 



exists, the synchronized state is globally asymptotically stable. Lower and upper bounds of Kl and Kjj 
are given by 

Kl> \ , and e:^<L£^, (17) 



5 = sup<'a min^^, a G M++ [> € (0, Anax/2]; (18) 



where 

a = suD < a 

a;e[0,2a] X 

Au = u — li-'*1ni with uj* defined in /max — ^^^x-^+oo /(icjfl;' ^max and A2(La) are the maximum 
degree and the algebraic connectivity of the graph, respectively. 

Proof. See Appendix Rl ■ 

Corollary 1 Assume that /(•), in addition to a2, is asymptotically convex or concav^. Then: 

1. If f{-) is unbounded, Kl = Kjj = 0; 

2. If f{-) is bounded, i.e., fmax < oo, upper and lower bounds of Ki and Kjj are 

i^L> V W ' """^ f WT ■ (19) 

Jmaxt^max Jmax^2 l.-'-'A j 

Remark 1. Even though conditions (|17|) or (|19|) provide only a range of values for Kl and Kjj, they 
state an important property of the whole system: If we want the network to reach a global consensus 
(common estimate), it is sufficient to take K greater than the upper bound in (jl7p or (jl9p : conversely, 
if we do not want the network to reach a global consensus, we need to take K smaller than the lower 
bound in (fT7|) (or (fT9]) ). In the simple case of unbounded coupling function, i.e. fmax = oo, Corollary [1] 
proves that the critical coupling value is -ftT = 0, since Kl = Kjj = 0. 

Remark 2. The nonlinear coupling model ([T|) includes, as a particular case, the linear coupling scheme, 
corresponding to the choice f{x) = x, as this function respects condition a2. From (jl9p . we realize that 
if the function f{x) is linear (and then unbounded), the critical value of K is zero. This means that a 
linearly coupled system always converges to a synchronized state, for any (positive) value of K (Corollary 
1). We may then ask ourselves whether there is any advantage in using a nonlinear as opposed to a linear, 
system. Indeed, the possibility to switch the system behavior from a system always converging to a global 
consensus to a system that cannot reach a global consensus, by acting on a single parameter, K, is a 
potential advantage that can be usefully exploited, as will be shown in the next section, to perform some 
kind of spatial clustering. This variety of behaviors is in fact a unique capability offered by nonlinear, as 



^The maximum of /(■) is defined on tfie extended real numbers, i.e., on R = R U { — oo, +oo}. 

function / : R R is said to be asymptotically convex or concave, if 33^ £ R : sign(/"(a;)) — sign{f" {x)),Wx > x, 

where sign(x) denotes the sign of x, and f"{x) is the second derivative of f{x) with respect to x. Observe that the above 

condition just avoids that /(•) could change its concavity infinitely often. Thus, it does not represent a strong restriction in 

the choice of the function /(■)■ 
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opposed to linear, systems. 

Remark 3. As a by-product of the proof of Theorem 1, in the particular case of c = 1, and under the 
conditions of Theorem 1, the dynamical system ([3]) approaches the synchronized state with a rate that 
is locally proportional to KX2(La)- The convergence rate is indeed a crucial parameter. In fact, from 
the point of view of the energy required to reach a common decision, the proposed system has several 
advantages with respect to centralized systems, as it is totally decentralized, but it has to pay these 
advantages with the energy wasted in the iterations necessary to reach the estimate. Clearly, the higher 
is the convergence rate, the lower is this waste of energy. The previous considerations suggest that, to 
increase the rate, we can increase K or change the network topology in order to increase A2(La), by 
increasing the network degree, for example. 

Remark 4. The wireless channel is typically affected by fading. Hence, it is important to analyze 
the proposed scheme when the coefficients aij are random variables. From the conditions required from 
Theorem 1, we realize that, provided that the nonzero fading coefficients are positive (this requirement 
has an impact on the kind of detector to be used) and the network is connected, the system maintains 
its capability to achieve a global consensus that is not affected by the values of aij. The only impact of 
the randomness of the coefficients aij is on the convergence rate, as A2(La) depends on them. 

The properties established by Theorem [H for scalar observation systems, can be generalized to the 
vector system ([5]), as follows. 

Theorem 2 Given system assume that conditions of Theorem{l\ are satisfied. Then, there exist two 
unique critical values of K, denoted by Kl and Kfj, with < Kl < Kjj, such that the synchronized state 
exists for all K > Kjj, and it does not for all K < K^. Furthermore, if it exists, the synchronized state 
is also globally asymptotically stable. Upper and lower bounds of Kl and Kjj are 

Kl > — — —, and Ku < rrr^: (20) 

/max"max g^2{i-'A) 

where g is defined in and Acj = [Aa;|", . . . , Alj^^] = Y'Y)q{u) — Itv ^ '^^), with u*^ given in (| i jp . 
Proof. See Appendix [Bl ■ 

Corollary 2 Assume that, /(•), in addition to a2, is asymptotically convex or concave. Then: 

1. If f{ ) is unbounded, Kl = Kjj = 0; 

2. If f{-) is bounded, upper and lower bounds of Kl and Kjj are 

maxi||Aa)i|| maxj ||AcJi||2 
Kl > — — —, and Ku < ' , ' • (21) 

5 Performance 

In this section we illustrate some properties of the self-synchronizing network proposed before. 



5.1 ML estimation in the presence of noise 

Different sources of noise affect the system: the observation noise, represented by the vector of random 
variables Wi in ([H]), and the system or couphng noise, represented by the stochastic process Vi{t) in ([1]) 
(or ([4])). These sources of noise affect system performance in a different way, as we show next. 

5.1.1 Observation Noise 

In Section 3, we showed how to choose the network parameters to guarantee the convergence of each 
node to the globally optimum ML estimate. However, in the case of unbounded noise, the convergence 
cannot be guaranteed with probability one. In fact, once we have chosen a coefficient K, there is always 
a non null probability that Kl > K, event that prevents the possibility of global synchronization. For 
any chosen K, using the upper bound in (llOp . the probability of the non-synchronization event can be 
upper bounded as follows: 

Pns := V{Kl >K}< V{Ku >K}< v!^\\D,Au\\2 > y/^axA2(L,)| . (22) 

In the case of Gaussian observation noise, Acj is a vector of i.i.d. zero mean Gaussian random variables 
and thus ||Z?cA^||2 is a x random variable, with degrees of freedom. Hence, denoting with Dz{z) 
the cumulative distribution function of Z := ||Z)cAa;||2, the probability that the network does not 
synchronize, for a given choice of K, is upper bounded as 

(K \ 

Pns < 1 — Dz ( ^ fmax^2{LA) 1 . 

This probability can be made arbitrarily small by using high values of K. 

As an example of vector ML estimation, in Fig. [2^), we show the derivatives 9i(t) as a function of time, 
for a network composed of = 16 sensors. The nonlinear coupling function /(x) in this case is bounded 
and equal to f{x) = tanh(x). Each node has degree 4. The observation model is the linear vector model 
(jl4p . with L = 3, M = 6; the matrices Aj are composed of i.i.d. Gaussian random variables of zero 
mean and unit variance. The common unknown vector is ^ = [1,2,3]"^. The observation noise is white 
Gaussian with zero mean and unit variance. The dashed line and the circles represent the global ML 
estimates achievable with an ideal control node that receives the observations and all mixing matrices 
Ai with no errors. We can see that, in spite of the low SNR, all the nodes converge to the ML estimate, 
as predicted. To get a global performance assessment, not conditioned to the specific realizations of 
the matrices Ai, in Fig. O we report the variances obtained in the same setting as in FiglJl averaging 
over 100 independent realizations of the mixing matrices Ai, the sensor initializations 9{0) and the noise 
values, as a function of the number of sensors A^. The network is a regular network, where all nodes 
have degree 4, for all values of A^. FiglS^) refers to the choice f{x) = tanh(x), whereas FigiSb) refers 




(dashed lines plus circles): a) observation noise only; b) observation plus coupling noise. 
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Figure 3: Estimation variance as a function of the number of sensors: a) f{x) = tanh(x); b) f{x) = sin(a;). 

to the choice f{x) = sin(x). The three sets of marks, in each curve, represent the variances obtained 
with the decentrahzed ML estimator, whereas the solid hnes refer to the centralized ML estimator. The 
convergence time is fixed to one second and it is the same for all N. It is interesting to observe that: 
1) even though sin(x) is not monotonic, and thus it does not satisfy assumption a2 of Theorem 1, the 
corresponding system behaves as the system with the monotonically increasing function tanh(x); 2) the 
decentralized method has practically the same performance as the centralized one; 3) even though the 
coupling is only local and it does not change with N, the variance decays as as the optimal ML 

estimator - this confirms the scalability of the proposed approach. 

5.1.2 Coupling Noise 

We focus now on the effect of coupling noise on the final estimate. Let us start with the effect of coupling 
noise on the conventional average consensus algorithms [HI [131 [HI [13 IS] • Without loss of generality, we 



consider as an example of average consensus algorithm, the discrete-time version of the linearly coupled 
dynamic system of [11^ \29^ 



where x[n\ = [xi[n\, - ■ ■ , a;7v[f^]]"^, with Xi[n] denoting the scalar state of i-th sensor at step n, and v[n] 
is the noise vector at step n and Wij is the weight associated by node i to the signal received from node 
j (Wij / if ttij 7^ 0, i.e., if nodes i and j are connected). If we pre-multiply (f23]) by l'^ and divide by 
A^, we get: 



with x[n\ = (l/N) ^^^^ This shows that the running average x[n] undergoes a random walk, thus 

implying that its variance increases linearly with the time index. This behavior was already observed 
in [29], whose authors realized that the average consensus achieved through (j23p does not converge 
in any statistical sense (except in the mean). Specifically, in j29] it was shown that, with average 
consensus algorithms, as given by (j23p . what converges to a constant value is the variance of the deviations 
Zi[n] := Xi[n] —x[n\. To recover from this problem, in [29] it was proposed a very elegant way to minimize 
the sum of the variances of Zi [n] , as a solution of a convex optimization problem, but the running average 
is still a process with variance increasing with time. 

Let us consider now the scalar system ([T]), but similar results can be obtained for the vectorial case 
Q . The study of the stability of the dynamical system ([T|) , in the case of nonlinear coupling and in the 
presence of noise (jH) is indeed a difficult problem and it goes beyond the scope of the present paper. 
Nevertheless, if we limit ourselves to the linear coupling case, i.e. f{x) = x, to make a comparison with 
the average consensus algorithm, as given in (j24p . we may still derive some basic properties. In fact, in 
the linear coupling case, exploiting the superposition principle, each dynamical system will converge to 
a random process having an average value, equal to uj*, as given in ([9|), corresponding to the solution of 
([T|), with uJi ^ and Vi{t) = 0, plus a fluctuating term having zero mean and variance cr^, corresponding 
to the solution of ([T|), with uji = and Vi{t) ^ 0. In other words, the effect of the coupling noise is to add 
a noise with constant variance, rather than of increasing variance, on the final estimate. This happens 
simply because the estimate is associated to the derivative of the state, rather than on the state itself. 

As a numerical example, in Fig. [2)3), we report a curve referring to the same settings as in Fig. [2^), 

except that now there is an additive white Gaussian process Vi{t) on each observed component, of zero 

mean and variance = 0.1. We can see that, also in this case, the state derivatives of all the nodes 

converge to values centered around the globally optimum ML estimates. 

^This model could also be applied, as the discrete-time counterpart of all the works considering average consensus through 
linear coupling. Furthermore, in the case of nonlinear coupling [13], (|23|l can be seen as an approximate version, valid when 
the states are close to each other and the additive noise is small. 



x[n] = Wx[n — 1] + v[n\,n = 1, 2, . . . 



(23) 




(24) 
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5.2 Effect of network topology 



From (jl9p . it is evident that the synchronization properties depend on the graph topology through the 
graph algebraic connectivity X2{La). This means that, for a given K and a given number of nodes, 
different topologies give rise to different behaviors. The easiest case to analyze is that of regular graphs, 
where all the node have the same degreqf|. The algebraic connectivity of an unweighted regular graph of 
degree d, having a ring topology where each node is coupled with only its neighbors, is 

A2 = 4^sm2(7ri/iV)p, ^-—^ (25) 

i=l 

where the last approximation is valid for d <^ N. In the more general case of non necessarily regular 
graph, the algebraic connectivity is not known in closed form, but it can be lower bounded as follows [33] 

A2>2(l-cos(^))<5-^5, (26) 



where 6 is the minimum degree and the approximation in (j26p is valid for ^ 1. Hence, in both cases, 
for a given d or 6, if the network size N increases, A2 decreases as l/N"^ . From (|19p . this means that to 
guarantee the self-synchronization, K must increase as A^^. 

Conversely, small world or scale-free random graphs exhibit a different behavior. Small worlds graphs 
exhibit, in fact, for any given degree, larger algebraic connectivity than regular graphs. As far as scale-free 
graphs are concerned, they are built starting from an initial number of nodes, say mg, and then adding 
new nodes according to an iterated procedure of growth and preferential attachment [36j. Denoting with 
^2("T'0 5 ^) the mean value of the second smallest eigenvalue of the Laplacian of a network composed of A^ 
nodes, averaged over the graph realizations, it was shown in |35j that the limit of A2(?tio, A^) for A^ going 
to infinity is constant. This proves that, if the network is built according to a scale-free topology, our 
strategy respects the scalability and fault tolerance properties, since, as soon as we choose K larger than 
the upper bound in (|19p . for sufficiently large A^, we are guaranteed that further addition or removal of 
a few nodes do not affect the global synchronization, and then estimation, capabilities of the network. 

5.3 Synchrony vs. desynchrony 

In Section 4, we showed that if the coupling is nonlinear and K is smaller than a critical value, the 
system does not converge. In this section, we show an example of the system behavior when we choose 
the coupling coefficient K in order to avoid the possibility to achieve global consensus. We considered 
a network of 40 x 40 sensors uniformly spaced over a regular planar grid. The initial measurements of 
the overall grid are reported in Fig. H] a) . Each sensor is initialized with its (noisy) observation and 
then it evolves according to ([1]). Each node is coupled with its neighbors and the maximum node degree 



The degree of a node is the number of its neighbors, i.e. the number of nodes Unked with that node. 



1 K 



is 12. A snapshot of the system state, after 1 second, is reported in Fig. [Hb). Comparing Figs. H] 
a) and b), we can see that the network is operating a sort of spatial clustering, even though the nodes 
are keeping evolving in time. The possibility of segmenting the observed field through a population of 




Figure 4: Intensity of the observed field: a) observed field; b) smoothed field. 



coupled dynamical systems is still an open research topic that we are currently investigating. 

6 Conclusion 

In this work we have shown that, if a sensor network observes a common event, a network of nonlinear ly 
coupled first-order dynamical systems can be used to achieve a globally optimum ML estimate, without 
the need to send any data to any fusion center. We have shown that the conditions guaranteeing the 
global asymptotic stability of the ML estimate, seen as the self-synchronization state of the whole sys- 
tem, depend on the coupling strength K and on the network topology through the algebraic connectivity 
A2. With respect to common average consensus techniques, based on the convergence of the state, the 
approach proposed in this work presents a stronger resilience against additive coupling noise. In general, 
the major advantages of the proposed strategy, as other average consensus or gossip algorithms, are the 
simplicity of each node, the scalability of the approach and the absence of congestion problems. These 
advantages are paid by the fact that the solution is achieved through an iterative procedure, whose 
convergence rate is proportional to the product K X2- In general, the total energy spent to achieve the 
final estimate, within a certain accuracy, is proportional to the product between the time necessary to 
reach the desired estimate (within prescribed error) and the power transmitted by each node. Since the 
iterative procedure involves only local exchange of information (provided that the overall connectivity is 
guaranteed), the transmit power of each node may be low. The time necessary to achieve the estimate 
is inversely proportional to K\2{L). To increase X2{L) it is necessary to increase the node degree, but 
this entails an increase of the transmit power of each node. Hence, we can foresee a sort of optimal 

1 a 



transmit energy, at each node, as a trade-off between tlie contrasting needs mentioned above. This is 
indeed an interesting research direction that we are currently investigating. For large scale networks, 
it is also important to study the effect of propagation delays. In [3^ we studied this problem and our 
results show that, for small delays, the system still converges, but the estimate becomes biased by an 
amount depending on the delays. Other interesting extensions include the effect of coupling noise for the 
general nonlinearly coupled system, the effect of having directed graphs, to model a system with different 
transmit power on each node, the effect of random coupling coefficients, to model channel fading effects, 
and the impact of time-varying topologies. 



Appendix 

A Proof of Theorem [1] 

We first introduce the following intermediate results, that will be used in the proof. 

Lemma 1 Given an oriented weighted graph ^ with N nodes, and positive numbers {wi}i associated to 
the edges, let — BD^B"^ be the (weighted) Laplacian ofW, where B is the N x E incidence matrix, 
Dw — diag (w) is the E x E diagonal matrix whose diagonal entries are the edge-weights Wi. Let Lw 
denote the generalized inverse of L^ f38^ . If the graph 5^ is connected, then lJv : 1^++ ^ is a 
continuous function in . 

Proof. Given w G consider the eigen-decomposition of L^ 

Lw = UwAwU^, (27) 
with eigenvalues arranged in nondecreasing order. The generalized inverse of Lw in ()27p is given by 

Lt = UwA;;^U^, (28) 

where Aw is the r x r diagonal matrix, containing the last r positive eigenvalues of Lw, and Uw the 
N X r matrix of the corresponding eigenvectors. Since the graph ^ is assumed to be connected, we have 
r = N — 1, i.e. |3H Lemma 13.1.1] 

rank(Aw) = N -1, Vw G M'^^. (29) 

The generalized inverse Lw in ([25]) is continuous in if it is continuous at w, for any fixed 

w G Rf^. Given w G Lw is continuous at w if, for any sequence {wj} of positive vectors 



Wj that converges to w, the corresponding sequence of generahzed inverses {Lw^jjgN converges to lJ^^ 
Stated in mathematical terms, we need to prove that 

V{w,}^.ep,-w {L«,^.},-eN-Lt, (30) 

where Lw^. denotes the generalized inverse of the weighted Laplacian L^j associated to the weights- vector 
Wj. Since is a continuous function of w, (jSOj) is equivalent to the following 

V{Lwj}jeN^Lw =^ {L^^.jjeN ^ Lw- (31) 

We prove now that (I3ip is satisfied, provided that ()29|) holds true (i.e., the graph is connected). 
To this end, we use the following necessary and sufficient condition for the continuity of the generalized 
(Drazin) inverse [38] Definition 7.2.3] 1^'^ . 

Theorem 3 ([38, Theorem 10.7.1][39| Theorem 2]) Let {Lw^ljeN Lw, with and given by 
[27\ ). Then, {Lw^ jjeN if and only if 



3jo > : rank(Aw,) = rank(Aw), Vj > jo- (32) 

If the graph is connected, from ([29]) it follows that rank(Awj) = rank(Aw), for all w,Wj- G ^++- 

Thus, given w G for any convergent sequence {LwjjjgN — > in (fSTI) . condition ([32]) is satisfied, 

and hence {Lw^ jjeN — > L^- This proves the continuity of at w G Since condition ([3T]) is satisfied 

for any w G Lw is continuous in all ■ 

Lemma 2 ( [40] Theorem 4.14]) Let C be a closed, convex subset of a normed linear space. Then, every 



compac 
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, continuous map F : C ^ C admits at least one fixed point. 



The proof of Theorem [Tj is organized as follows. We introduce, first, a proper transformation of the 
original system ([3]), so that the existence and the global asymptotic stability (according to Definition 
d]) of the synchronized state ([9]) can be recasted in the classical study of existence and the asymptotic 
stability of the equilibria of the transformed system (see, e.g., |41]l42j). Then, using standard fixed point 
arguments, we prove that, under al-a3, an equilibrium for the transformed system exists, provided that 
K > Kij^ and cannot exists \i K < K^. Finally, we show, introducing a valid Lyapunov function, that, 
if an equilibrium exists, it is also asymptotically stable. 



sequence {Aj}^gj, of matrices in R"^" is said to converge to A gR"^", if for every real number e > 0, there exists 

an index jo(e) such that if j > jo(e), then ||Aj — A|| < e, where |[-|| is some matrix norm on R"^". We denote a convergent 

sequence {Ajj^gj^ to A as {Ajj^gj^ A. It is worth observing that {Ajj^gjj A if and only if the entries of Aj converge 

to the corresponding entries of A. 

^"For the sake of simplicity, we adapt [38| Theorem 10.7.1] to our notation. Observe that the Drazin inverse, as defined 

in [381 Definition 7.2.3], corresponds to the generalized inverse given in (|28|l . 

^^The map / : C t— > C is called compact if /(C) is contained in a compact subset of C. 

t Q 



A.l Existence of the synchronized state 

Let us assume that conditions al-a3 are satisfied and consider the following change of variables 



^i{t) = 9^{t)-uj^t, i = l,2,...,N, (33) 

with Lo* defined in ([9]). The original system ([3]) can be equivalently rewritten as 

= Ao; - KD-^BDa/ [B'^*(t)] (34) 

= Au- KD~^BDAD*B^*(t) (35) 

= Ao; - ^LA,**(t), (36) 



where *(t) = [^'i(t), . . . , ^'7v(t)]^ , with *(0) = 0(0), Au =uj — u>*l, and La,* is the weighted Laplacian 



of the graph, with diagonal weights-matrix DaD^ (that depends on and [D^] •• is given by 

[Di-k= Lr.J' >0^ V*GM^, i = l,...,E, (37) 



where the positivity of [D^] •• > comes from alF^. 

According to Definition 1, the synchronized state of ^ exists if and only if the dynamical system in 
(j36p admits an equilibrium, or equivalently, if there exists a solution for the following system of nonlinear 
equations 

La,** =-^DeAc^. (38) 



Given (j38p . we prove that there exist two critical non-negativa I values of K, denoted by Kl and 
Kjj, such that for all K > Kjj the system (j38p is feasible and for all K < Kl the solution of (j38p 



disappears. To this end, it is sufficient to provide a lower bound of Kl, and an upper bound 

K^^^ of Ku, such that for all K > K^"^^ the system ()38p admits a solution, and for all K < K^"'"^ does 
not. Given these and K^'^\ there must exist a unique Kl and Ku as defined above such that 

^(low) <Kl<Ku< kJ^p^ 

It is worth observing that, in the case of the linear function /(x) = x, we have [D*]-- = 1 and 
La,* = La, V'J' € M^. Hence, the system (|38p becomes linear and, since l^DcAci; = and LA,*liv = 
Oat, it admits, for any K ^ 0, oo'^ solutions, given by ^* = L^D^Aa;/K + span{l7v}, where L^ is the 
generalized inverse of the weighted Laplacian La [38]. Thus, in the case of f{x) = x, we have Kjj = 
(by definition) and Kl = (since for K < the equilibrium points of (j36p are not stable, as shown in 
Appendix |A2]). 



^^Note, from item 3, right after (1), that lim^i-^o f{x)/x = /(O) = 1 
^We focus only on nonnegative K, since the potentia' 

are not stable for ([3]), as we will show in Appendix lA. 21 



^^We focus only on nonnegative K, since the potential solutions of the system (|38p corresponding to negative values of K 



in 



Conversely, in the case of nonlinear bounded functions /(•), i.e. lim2,.^+oo f{x) = /max < +00, is 
lower bounded by a positive quantity. This lower bound corresponds to the value of K below which the 
system (138]) is surely infeasible, i.e. (see (p4|l ) 



K ||BDa/ [B^*] IL < llDcAcc^ll^ , (39) 
where \\-\\^ denotes the infinity norm of a vector. A more stringent condition than (j39p is 

i^C?max/max < ||DcAcj||^ , (40) 

n A ^ 

where we used the inequalitjo HBDaHo^^ < dmax, with c?max — max, [BBAlij ■ From (gOD it follows 

j=i 

that, if the system (j38p admits a solution, then it must be 

K > 4 ^^^J^, (41) 

Omax /max 

which provides the lower bound in (fT7|l . Observe that, for unbounded nonlinear functions (i.e. /max = 
+oo), the lower bound ([^T|) disappears. 

We consider now a generic (bounded or unbounded) nonlinear function /(•), and provide sufficient 
conditions on K for the system ()38p to be feasible. A solution for (138p exists, if and only if the following 
mapping admits at least one fixed point in 

* =^^DeAa;, (42) 
where ^ is the generalized inverse of the weighted Laplacian La,*, given by 

-'-'A,* = Ua,*Aa^*Ua,-*, (43) 

with Aa,* the {N — 1) x (N — 1) diagonal matrix, containing the last N — 1 positive eigenvalues 
of La,* (assumed to be arranged in nondecreasing order), and Ua,* the x (A^ — 1) matrix of the 
corresponding eigenvectors. To prove the existence of at least one fixed-point for (I42p it is sufficient to 
show that (|42p admits a solution in some compact, convex set of M^. We choose this set, without loss of 
generality (w.l.o.g.) as = {* G : ||*||2 < a} , where a is any positive number Observe that 

Ba is compact and convex for all a G M-f+, and that, by Lemma [H the mapping L^^ Au/K in (I42p is 
continuous on M.^ (because of the positivity of [D*]^,- > 0, G M^). Hence, according to Lemma O a 
fixed-point for (j42p exists in Ba, if L^ ^D^Acj/AT is a compact map on Ba, for some given a € This 
is guaranteed if, for any given a G K is chosen so that ||L^ ^DcAa;||2 < K a, which corresponds 

*° « 

K > ± = " / V* G Ba, and a G (44) 

a a \2 (La,*) 



The matrix norm induced by the vector infinity norm is the maximum among the absolute values of the row sums [381 
Proposition 10.2.2]. 



on 



where ||L^ ^||2 is the spectral norm of and A2 (La,*) is the second smallest eigenvalue of La,*- 

To remove the dependence of A2 (La,*) on we consider the more stringent (sufficient) condition 

< . ' , X r <K, a G M++, (45) 



«A2(La,*) amin^-gCa A2 (La,* 
where we used the inequality min^gCa -^2 (La,*) < miii*eBa A2 (La,*) , with Ca = [— a,a]^ 5 B, 
Theorem 1.16]. In (j45p . the minimum of A2 (La,*) on Ca can be lower bounded as follows. For all 
* > * > 0, since L^-^ - L^ ^ ^ 0, we have 144, Problem 4.3.14] 

Ai(L^^^) > Ai(L^^^), V*>*>0^, z = l,...,iV, (46) 

where Aj(L^-^) and Aj(L^^) denote the i-ih. eigenvalue of L^-^ and L^^, respectively, arranged 
according to the same order. From (06]), it follows that the lower bound of the minimum of A2 (La,*) 
on Ca occurs for the minimum of the weights [D^]-- , defined in (|37p . Since [B-^^f] . G [—2a, 2a] for any 
^ £Ca, we have 

A2(La,*)>A2(La) min M = A2 (La) min V* G C„, (47) 

xe[-2a,2a] X x6[0,2a] X 

where the last equality follows from the fact that, because of al, f{x)/x is even. Using (j45p and ()47p . 
we obtain the following bound for K : 

A- > A» A „,R,„ (48) 

amma,g[o,2a] 

For any given a G K{a) defined in ()48p represents the smallest value of K, for which a solution 

of (j42p is guaranteed to exists in Ba- Increasing a G we enlarge the region Ba where the solution 

may fall. Since any K{a) in (jl8]) . with a G ^++, is a valid upper bound for Ku, the lowest upper bound 
of Kjj is obtained taking the greatest lower bound (i.e. the infimum) of the set 

JC = {K{a) : a G M++} , (49) 

with K{a) given by ()48p . In fact, by the approximation property of the infimum [431 Theorem 1.14], for 
every K > inf/C, there always exists some K{a) G /C such that K > K{a) > inf/C, which, by (j48p . is 
sufficient to guarantee a solution of (02]) . Since /C in ([19|) is bounded from below and non-empty, inf/C 
always exists [43, Axiom 10] and is given by 

inf^=^^, (50) 
5A2 (La) 



whert 
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I /(^) 1 
= sup < a min , a G M++ > . (51) 

x€[0,2a\ X ' 



^^The matrix norm consistent with the Euclidean vector norm is the spectral norm [38i Proposition 10.2.4], defined as 

the largest singular value of the matrix. 

^^Note that the sup in (|51[) is defined on the extended real numbers, i.e. on R = R U {— cxd, +00} . 

01 



Since min^g[o,2a] f{^)/^ ^ /(2a)/(2a), Va € expression ([5T]) can be upper bounded by [43} Theorem 
1.16] 

/ • ■^(^) ^ TO 1 ^ / •^(^") ^ TO 1 -^max 

5 = sup < a mm , a G > < sup < a— — , a € M++ > = , (52) 



xe[0,2a] X J 1^ 2a 

which provides the following lower bound for inf /C in (j50p : 

inf/C> " ; (53) 

/max -^2 V-LiAj 

This complete the proof of (fT7|) . 

Observe that the lower bound in (I53p may be reached or not, depending on the particular function 
f{x). Without additional properties on /(•), we have no guarantee about the achievability. 

We prove now that a sufficient condition for inf /C to satisfy (j53p with the equality is that /(•) be 
asymptotically convex or concave. Stated in mathematical terms, a function /(■) is asymptotically convex 
or concave if 

a4) : 3 X G M : sign (/"(x)) = sign (/"(x)) , Vx > x, (54) 

where sign(x) denotes the sign of x, and /"(x) is the second derivative of /(x) with respect to x. 

Under al-a4, the function /(x)/x is continuous on R and it is quasi-convex or quasi-concave on 
[x, +cxd), where x is defined in ()54p . In fact, because of a4, /(x) is convex (or concave) on [x, +oo), and 
thus the set {x € [x, +oo) : /(x) — ax < 0} (or {x G [x, +oo) : /(x) — ax > 0} ), with a G M, is convex, 
because it is the sublevel set (or the superlevel set) of the convex (or concave) function /(x) — ax; which 
corresponds to the definition of quasi-convexity (or quasi-concavity) of /(x)/x on [x, +oo) [45j. Hence, 
one of the two following statements must hold for /(x)/x on [x, +oo) [45., Section 3.4.2 

3 X € [x, +oo) : Vx € [x, +oo), /(x)/x is nondecr easing, (55) 

if and only if /(x)/x is quasi-convex on [x, +oo), or 

3 X G [x, +oo) : Vx € [x, +oo), /(x)/x is nonincreasing, (56) 

if and only /(x)/x is quasi-concave on [x, +cx3). Using (f55l) and (f56l) . we obtain the following result for 
the minimum of /(x)/x on [0, 2a 



If /(x)/x is quasi-convex on [x,+oo), then (see ([55]) ) 

_ /(x) f(w) 
3 a S [x,+(X)) : min = , Va E [a/2,+oo) and some w G [0, a]; (57) 

a;G[0,2a] X W 



^'^The class of functions f{x)/x that are nonincreasing (or nondecreasing) on the whole interval [a;, +oo] is tacitely treated 

as special case of (|55p and (|56|l . as shown after (|57p and (|58|) . 

'^^We assumed, w.l.o.g., that a > 0, since, if a < 0, conditions (|55p and H56[) are satisfied for any j: > 0. 



oo 



If f{x)/x is quasi-concave, either (j57p or the fohowing condition hold true (see ([561 



3aG[x,+cx3): mm = — ^, Va G [a/2, +00). 



(58) 



x£[0,2a] X 2a 

Observe that, in the case where f{x)/x (satisfying (|54p ) is nondecreasing (nonincreasing) on the 
whole interval [x, +00), condition ([57]) (condition ([58]) ) still holds true. For a bounded function /(•), 
lim^^+00 f{x)/x = 0, and thus (only) condition ([58|) is satisfied. In the case of unbounded /(•), instead, 
either ([57|) or ([58[) may be met. 

Using ([57]) and ([58]) . we show that, under a4, inf defined in ([50|) achieves the lower bound in ([53[) . 
To this end, consider the following subset of IC 



ICa = {K{a) : a € [a, +00)} C IC, 



(59) 



where i^'(a) is given in (|48p and a is implicitly defined in (|57p . if /(x) is asymptotically convex, or in 



58]) . if f{x) is asymptotically concave 
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K{a) = < 



Given (1571) and (l58l). i^(a) G /C;? can be written as 



Va G [a/2, +00). 



if (E3) holds, 

a A2 (La) 



2||DeAa;||2 if ([58]) holds, 



(60) 



t /(2a)A2(LA)' 

Whether ([57l) or ([58l) is satisfied, -ftr(a) in ([60|l is a continuous and decreasing function on [a/2, +00). It 
follows that 



inf /C^r = lim K(a) 



a — ' + 00 

a>a 



0, 



2 ||DcAw||2 

. /maxA2 (La) 



if /(•) is unbounded. 



if /(•) is bounded. 



(61) 



Using ([53]), ([59]) and ([H]) we obtairP 

2 llDpAo; 



lo 2||D,-Alj|L 
^ < inf/C <inf/C;r - " 



(62) 



/maxA2 (La) " /maxA2 (La) ' 

where the second inequality in ([62]) comes out from /C5 C /C. From ([62]) it follows that inf /C =inf/C2, 
which proves the upper bound in (|19p . 



A. 2 Global stability of the synchronized state 

Assume now that, in addition to conditions al and a3, K > Kjj > 0, so that system ([3]) may synchronize. 
We prove that the synchronized state of the system ([3]), whose existence is guaranteed hy K > Ku, is 
globally asymptotically stable (according to Definition [1]). 



^^With a slight abuse of notation, we use the same symbol a for both of the conditions (|57|) and (|58p . 

■^"We use a unified expression for inf/Ca, for both bounded and unbounded functions /(■), with the convention that, if 

/max = +00, then inf /Ca = 0. 
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To this end, it is sufficient to consider system p4p and sliow tliat the state of (|34p converges to an 
equihbrium, from any set of initial conditions. Left-multiplying (I34|) by l^Dc and using ([2]), we obtain 

l?^Dc*(t) = l?^DcAa; - K1?^BDa/ (B^*(t)) =0 l?^Dc*(t) = l^Dc*(0), Vt > 0. (63) 

In words, the weighted sum l^Dc'J'(t) is an invariant for system (j34p . The invariance of l^Dc'J'(t) 
allows the following decomposition of the state vector ^{t) in ([3i 



with 

l?;DeA(t) = 0. (65) 
Thus, we can study the evolution of system (I36p by studying the dynamics of the following system 

A{t) =Au- ^BDa/ (B^A(t)) , (66) 

with A(t) satisfying the constraint (j65p . 

We focus now on the equilibria of (|65p and show that system (|66p admits a unique equilibrium and that 
such an equilibrium is globally asymptotically stable; this proves also the globally asymptotic stability 
of the synchronized state of ([3]). 

Since c^Iat ^ 0, the condition K > Kjj guarantees the existence of an equilibrium for system ()66p 
(c.f. Appendix lA.ip . Moreover, because of ()65p . all the equilibria of (I66p are isolateJ^. In fact, the 
Jacobian of BDa/ (B-^A(i)) is given by BDa diag(/(B-^ A))B^, which is positive definite for all the 
vectors A that are equilibria of (|66p . Denoting by A* = [A*, . . . , A^]^ one of the (isolated) equilibria 
of ([66]) . i.e. satisfying 

Au = KB-^BBaJ (B^A*) , (67) 

by translating the origin to the equilibrium A'^, we can make Oat to be an equilibrium of (I66p and writJ^ 

A{t) = -i^D-^BDA [/ (B^A(t) + B^A'^) - / (B^A'^)] . (68) 

To make explicit the dependence of ([68]) on the weighted Laplacian, we introduce the following 
function 

OA^M^) ^f{x + A)- A*) - / (A^^ - A^) , (69) 
that, because of al, satisfies the following properties 

( dfix) 



9a* A* (^) 



X 



dx 

> 0, X / 



a;=A*-A* 

3 » 



> 0, X = 0, VA* : A* ^ A*. (70) 



^^An equilibrium point is isolated if it has a surrounding neighborhood containing no other equilibria. 

^^With a slight abuse of notation, we use the same variables A, to denote also the system (|66|l . after the shift around A* 

OA 



Using ([69]) and introducing the diagonal matrix D^*^^, whose positive (see ([70]) ) diagonal entries are 
given by gA*,A*{^j{t) — Aj(t))/(Aj(t) — Aj(t)), indexed from 1 to E, and 5a*,a*(') is defined in (f69]l . 
system (I68p can be equivalently rewritten as 

A(t) = -KD-iBDADA*,AB^A(t) 

^ -i^D-iLA,A*,AA(t), (71) 

where a* a is the weighted Laplacian of the graph, with diagonal weigth-matrix DaDa* a- 

The system (|71|) admits the point O^r as the unique equilibrium (which also guarantees the uniqueness 

of A* for (j66p ) and such an equilibrium is globally asymptotically stable, as we argue next. 
From al and the properties of the weighed Laplacian of a connected graph, we have that 

La,a*,aA = Ojv ^ a G span{ljv}, VA*. (72) 

From (|65]) and (f72|l . it follow that all the equilibria A of system ([71]) are given by 

A G {span{lAf} n span{c_L}} = {On}, (73) 

where spanjcx} = {x E M^: l^DcX = 0}. Hence, the unique equilibrium of ()7ip is the vector A = Oat. 

After having showed the uniqueness of the equilibrium, we can prove its global asymptotic stability. 
To this end, consider the following candidate positive definitj^ Lyapunov function 

V{A{t)) = \ II Dy^A(t) f . (74) 

The function y(A(t)) in ([71]) is non-increasing along trajectories of ([7T]) . since 

y(A(t)) = A^(t)DeA(i) 

= -i^A^(t)LA,A*,AwA(t) < 0, (75) 

where the last inequality follows from LA,A*,A(t) ^ OnxN, VA*, A(t)G M", and equality in ([7^ is reached 
if and only if La a* a A = Oat. Since, by definition, A must satisfy also the constraint l^DcA = (see 
([65]) ). the function ^(A) = if and only if ([73]) holds true, i.e. A = Oat. Hence, V{A) in ([H]) is a valid 
Lyapunov function for ([7ip and A = Oat is globally asymptotically stable for (|7ip [42] Theorem 5.24]. 
This guarantees that the state ^{t) of ([M]) converges to an equilibrium of ([M]) as t — > +cxd, for any set 
of initial conditions. 

Observe that, through the whole proof, we have always considered positive values of K. This comes 

from the fact that, all the potential equilibria of the system p4p (and thus of ([7ip ) corresponding to 

^^A continuous function V : K" R is called a positive definite function if V(0) = 0, and V^(x) > ci(|x|), Vx G R", where 
Q : R 1-^ R is some continuous, strictly increasing scalar function, with a(0) = 0, and a(p) h-> oo as p oo [421 Definition 
5.13]. 
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negative are instable. In fact, for negative the valid Lyapunov function y(A(t)) defined in ()74p . 
has first derivative V{A.[t)) along trajectories of (I7ip that is positive definite, which proves the instability 
of the equilibrium A = Oat [421 Theorem 5.29]. 

B Proof of Theorem [2] 

Since most of the proof of Theorem follows the same approach, already described in the Appendix El in 
the following we point out only the differences. 

B.l Existence of the synchronized state. 

Given the system ([5]), we introduce, as for (|34p . the following change of variables 

^^{t) = ei{t)-ult, i = l,2,...,Ar, (76) 

with cj^ defined in (jlip . and rewrite ^ as 

= - K-D^V^ilL ® BDA)f {{II ® B^)V^{t)) , (77) 

where *(t) = [*i(t), . . . , *Ar(t)]'^ , with *i(0) = 0^(0), and I\(X)=u - {In ® ^i)- To obtain ([TZD, we 
have used the following chain of equalities P(lAr ® u:*j^) = (cj^ (g) l^r), and {II ® B^){u*j^ (g) l^r) = 
diag(a;^) ® B'^In = OlexL- Given ([77]) . the synchronized state of ([5]) exists if and only if the following 
system of non-linear equations admits a solution 

PDqAu; = KilL BDa)/ {{II » B^)P*) ■ (78) 

We recast now the study of the existence of a solution for (j78p , to the study of the solution of a set of 
simpler sub-systems, similar to (j42p . To this end, we introduce the vectors ^ = P^' and Acj = PDq Au;, 
partitioned as * = [*^, . . . , *^]^, and Alj = [Acjf, ...,Aulf, with *i = . . . ,^^j!^f and 

= [Auj^^\ . . . , Au;^"*]"^ , i = 1, . . . , L, where ^j*'' denotes the i-th component of the j-th. sensor's 
state vector ^'j, and ^^^j^ is the i-th component of the vector Uj — cj^. Then, the system (j78p can be 
equivalently rewritten as 

L*i,A*i =;^^'^*> i = l,...,L, (79) 

where L^. ^ = BDaD^.B^ is the weighted Laplacian of the graph, with diagonal weights-matrix 
DaD^., and D^. is still given by (iST]) . by replacing ^' with ^'j. Thus, a solution of ([751) exists if and 
only if every equation in ()79p admits a fixed point in R^. But, for any given i, ()79p is equivalent to (I42p . 
if La.* and DcAo; are replaced with L^. ^ and Acjj, respectively. Hence, each equation in (j79p admits 
a fixed point provided that K > whereas a solution of (j79p cannot exist if < i^'^*^ A lower bound 



(i) (i) 

for Kj^ and an upper bound for Kjj are given by 



Af>|^.. a„d K<J><:^ (80) 



with g defined in (jlSp . 

Since the couphng among the equations in ()79p is given only by the presence of K, a fixed point for 
all the equations exists if and only if the critical values and Kfj in (I79p are chosen as the maximum 
of {K^ji and {Kjj^}i respectively, with and Kjj'^ defined in ([801). This proves (|20l) . 

The upper bound of i^^/ and lower bounds of Kl in (|21|) can be obtained following the same approach 
used in Appendix lA.il to prove p9]) . 

B.2 Global stability of the synchronized state 

We consider the system ()77p and show that, under al — a3 and K > Kjj, the state vector converges to 
an equilibrium, from any set of initial conditions, which proves the globally asymptotic stability of the 
synchronized state of ([5]). Since, similarly to ([63j) 

N ^ N N 

(l^®lL)DQ*(t) = J];Qi*,(t) = 0l ^ ^Qi*i(i) = J];Q,*i(0), Vt>0, (81) 

i=l i=l 1=1 

[ij^ Ij^)DQ^'(t) is an invariant and thus we can decompose ^{t) in (j77p according to 

*(t)^(ljVxJV®lL)DQ*(0) + A(t), (82) 

where 1 AT X Af is the iVxiV matrix of all ones, Dq = (in <S) (Sili Qi) DQ,andA(t) = [AJ (t), . . . , Al,{t)f 
satisfies the constraint 

TV 

(l?^0lL)DQA(t) = ^QiA,(t) = OL. (83) 

i=l 

Introducing (j82p in (j77p . we obtain the following dynamic for A(t): 

A(t) = Ac^ - iCDQip^(lL ® S£>a)/ ((i"L ® S^)PA(t)) , (84) 

where we used the chain of equalities {II B'^)I'{1nxn ® II) = {II ^ -B"^)(Il ® lAf)(l|i- <S5 II) = 
(/l ® B^1n){1% «) II) = OiijxLTV. 

Following the same approach used in Appendix IA.2I to obtain (|68p . we can translate the system ()84p 
around the isolated equilibrium A* given by 

Au = KDq1p^(/l S£)a)/ ((/l ® S^)PA*) , 

(85) 

(l^®Ii)DQA^ = OL, 

so that the vector Oln is an isolated equilibrium of the following translated systeiq^ 

A(t) = -i^DQip^(/L BDa) if {{II ® S^)P(A(t) + A*)) - / {{II S^)PA'^)) . (86) 



*We use the same variables to denote the translated system (|84p around A* 



We rewrite now the system (j86p in a more compact form, making explicit the dependence on the 
weighted Laplacian. To this end, we introduce the vectors A(t) = PA(t) and A* = PA*, parti- 
tioned as A(t) = [A'^{t), Alit)f, and A* = [Af , . . . , Aff, with A,{t) ^ [A?{t), A^^{t)f 
and A* = [A^^*'',... jA^*^]"^ , where A^j\t) and A*^*^ denote the z-th component of Aj{t) and A^, 
respectively. Then, the system ()86p can be equivalently rewritten as 

A(i) = -KD-ip^LA,A*,APA(t), (87) 

where A{t) satisfies the constraint ([83l) . and L a, A*, A — i^L ^ -B)Da,a*,a(-^l ^ -S'^)) with Da. a*. A — 
(Il'Si'Da) diag(D^* ^ , • • • , D^* ^ ) and D^* ^ is the Ex E diagonal positive matrix, whose diagonal 
entries are given by g^i,\k) ^*(k){A^p — A^^'')/{A^p — A^^) indexed from 1 to E, and g^^,(k) ^*(fc)(-) is the 
same function as defined in ()69p . with Aj and A* replaced by AJ and A^- , respectively. 

Using the same technique as in Appendix IA.2I one can prove that the vector Olat is the unique 
equilibrium of (jSTj) and it is globally asymptotically stable. A valid positive definite Lyapunov function 
for dHZl) is 

V{A) = 1/2 II D^/2A(t) ||2, (88) 

that can be seen to be non-increasing along trajectories of (j86p . and with zero derivative if and only 
A = Oln- This proves the globally asymptotically stability of the synchronized state of ([5]). 
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